Competing associations in six-species predator-prey 
models 



Gyorgy Szabo 

Research Institute for Technical Physics and Materials Science, P.O. Box 49, H-1525 
Budapest, Hungary 

Abstract. We study a set of six-species ecological models where each species has two 
predators and two prey. On a square lattice the time evolution is governed by iterated 
invasions between the neighboring predator-prey pairs chosen at random and by a site 
exchange with a probability Xg between the neutral pairs. These models involve the 
possibility of spontaneous formation of different defensive alliances whose members 
protect each other from the external invaders. The Monte Carlo simulations show a 
surprisingly rich variety of the stable spatial distributions of species and subsequent 
phase transitions when tuning the control parameter Xs- These very simple models are 
able to demonstrate that the competition between these associations influences their 
composition. Sometimes the dominant association is developed via a domain growth. 
In other cases larger and larger invasion processes precede the prevalence of one of 
the stable associations. Under some conditions the survival of all the species can be 
maintained by the cyclic dominance occurring between these associations. 



1. Introduction 

Ecological systems exhibit a large number of associations of certain species which coexist 
on the same territory. These associations can be considered as evolutionary stable states 
(or Nash equilibria) whose stability and competitiveness depend on the environmental 
conditions [H El El IH El IE] • Generally, the systematic analysis of the competition between 
these spatial states is made difficult by the complexity of the phenomena and by the 
large number of parameters in the suitable models. 

In order to have a deeper insight into these phenomena we consider now a family 
of toy models exemplifying some basic features and the richness of possible behaviors. 
In these predator-prey models the species, located on the sites of a lattice, invades one 
of the neighboring sites if it is occupied by its prey. The investigation of predator- 
prey systems on lattices has a long history reviewed recently in the papers [SI 13 Ej- 
In these systems the invasion processes can yield extinction of species and/or maintain 
some steady states in the spatial distribution if the corresponding food web (a directed 
graph) contains loops. In such a situation a self-organizing spatio-temporal pattern can 
be formed by the species dominating cyclically each other. Besides it, in the present 
model the neutral (neighboring) species are allowed to exchange their sites. This latter 
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process helps the formation of many well-mixed distributions of neutral species which 
can also protect each other against the rest of species. As a result these models possess 
many final states toward which such a system can evolve. In the spatial systems these 
states can coexist by forming domains if the system is started from a random initial state 
for sufficiently large sizes. The competition between these multi-species associations 
along the boundaries can affect their internal structure and the evolutionary process 
selects one of these states. Thus, the investigation of these simplified models can help 
us to clarify some mechanisms related to Darwinian selection and the emergence of 
structural complexity [HI characterizing pre-biotic evolution too |lUj . 

Several basic features are already well known for the simplest versions of the 
above models. As mentioned above, if the spatial predator-prey model on the square 
lattice contains three species with cyclic dominance (this model is also called as spatial 
rock-scissors-paper game) then the cyclic invasions maintain a self-organizing pattern 
with equal species concentrations jTH [T21 US] • Both the composition [T31 US] and the 
geometrical features jTHl Ej of this spatial pattern are affected by the independent 
variation of invasion rates meanwhile this structure provides a stability against some 
external invaders fHl Ej- Similar self-organizing pattern is found for those cyclic 
predator-prey (or voter) models where the number A^^ of species is limited, namely 
3 < A^s < 14 113. For even A^^, however, these systems are very sensitive to the 
independent variation of invasion rates jSOj- More complex behavior was observed for 
those systems where the food web consists of several cycles that can be distinguished 
from a topological point of view ^H]- In these multi-cycle ecological models the final 
(stationary) state is prevailed by a defensive alliance being the most stable cyclic 
subsystem. These defensive alliances are composed of those species that form a cycle 
in the food web thereby their spatial distribution is equivalent to those of one-cycle 
systems mentioned above. The privileged role of these associations is related to their 
most effective protection, provided by the topological features of the food web, against 
the external species. In their enhanced protection the role of spatial structure is crucial 
because under mean-field conditions these associations are not able to eliminate the 
external invaders. 

Very recently the investigation of a four-species cyclic predator-prey model has 
indicated the spontaneous formation of another type of defensive alliances [21]. In this 
latter model the lattice sites can be empty and the individuals are allowed to jump to 
one of the neighboring empty sites. As a result the odd (even) label species can form 
a well-mixed state ensuring the possibility of prompt counter-attack against the even 
(odd) label species and by this means these associations can defend their territories. 
In other words, the neutral species can also form defensive alliances in the presence of 
some local mixing. 

In the present paper we consider a set of models allowing the appearance of both 
types of above mentioned defensive alliances. In these six-species models the food web 
contains many cycles and the local mixing is provided by a site exchange between 
the neutral pairs. For all the models the probability Xg of site exchange will be the 
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only tunable parameter. Using Monte Carlo (MC) simulations it will be shown that, 
despite their simplicity, these models are able to account different stationary states and 
subsequent phase transitions when varying the value of Xg. The stationary states are 
composed of different species whose number varies from 2 to 6. Evidently, the time scales 
of the formation of these states are dependent on the microscopic mechanism (invasion 
or mixing). Sometimes equivalent states are formed via a nucleation mechanism and 
they can destroy each other (cyclically) when the boundaries of these growing domains 
meet. 

Apparently the combination of our previous models 1^ means a simple 
extension of the six-species models. This extension, however, exemplifies that a small 
variation of a sufficiently complex ecological model can yield surprisingly large increase 
in the types of behaviors and phase transitions whose clarification raises many other 
questions. At the same time, many aspects of these features seem to be valid for 
the sufficiently complex ecological systems (or evolutionary games) as discussed in the 
Conclusions. 

2. Spatial predator-prey models 

In all the present models the site i of a square lattice is occupied by an individual 
belonging to one of the six species and the corresponding distribution is given by a 
set of site variables (sj = 0, ... ,5). The predator-prey relation between the species is 
defined by a food web. On the corresponding directed graph the nodes (labelled from 
to 5) refer to the species and the directed edges point to the prey from its predator. 
Our analysis is focused on those systems where each species has two predators, two 
prey, and a neutral one (not linked). In this case each possible food web (satisfying 
the above conditions) is isomorphic to one of the four graphs shown in Fig. ^ and the 
corresponding model is named after their labels A, B, C, and D P^ . 

The evolution of the species distribution is governed by iterating nearest neighbor 
invasions and site exchanges. More precisely, starting from a random initial state we 
repeat the following steps: 1) two nearest-neighbor sites are chosen at random; 2) if 
these sites are occupied by a predator-prey pair then the prey is killed and an offspring 
of the predator occupies the site of prey; 3) for neutral pairs the species exchange 
their sites with a probability Xg. Evidently, nothing happens if the individuals belong 
to the same species. For example, in the model A the randomly chosen pair (1,2) 
transforms into (1,1) (invasion) and the pair (1,4) transform to (4,1) (site exchange) 
with a probability Xg, otherwise the neutral pair remains unchanged. Notice that the 
invasion rates between the predator-prey pairs are unique and the models have only one 
parameter (Xg) characterizing the strength of mixing between the neutral species. 

These models are investigated by a series of MC simulations performed on a square 
lattice with sites N = L x L under periodic boundary conditions at different values of 
Xg. The simulations are started from a random initial state. Within a time unit (MCS) 
each pair has a chance once on the average to perform an invasion (between predator 




Figure 1. Food webs defining the four classes of predator-prey relations where each 
species has two predators and two prey. The decorated lines connect the members of 
the defensive alliances providing mutual protection against the external invaders via 
cyclic invasions. 

and prey) or a site exchange with a probabihty Xs (between neutral species). During 
the simulations we have recorded the current values of species concentrations [pait)', 
a — 0, . . . , 5] as well as the probability of finding predator-prey [-Ppp(i)] and neutral 
pairs Pn{t) on two neighboring sites. The stationary states are characterized by the 
averages values of these quantities (denoted as p^, Ppp, and P„) which are determined 
by averaging over a long sampling time varied from 10^ to 10® MCS after a suitable 
thcrmalization time. Evidently, Pa = ^ and the quantity 1 — Ppp — Pn gives the 
probability of finding the same species on two neighboring sites. At the same time we 
have determined the fluctuation of these quantities deflned as Xa — ^[{Pa) ~ (Pa)^] 
where (...) refers to the averaging over the sampling time. To avoid undesired size 
effects the hnear size is varied from L — 500 to 2500. The larger sizes are used in the 
close vicinity of transition points where the fluctuations diverge. 

Now we discuss the general featiu^cs of the possible stationary states. For all these 
models there exist six homogeneous states denoted as ^^^-^ where Si = a (a = 0, . . . , 5) 
for each site i. These states remain unchanged and the cvohition of the species 
distribution is stopped whenever one of these (absorbing) states is reached. At the same 
time, the homogeneous state is unstable against the invasion of the suitable predators. 

The mixed states of neutral species are denoted as $03^ if Sj = or 3. Evidently 
there exist two additional neutral pairs for all the four models. In such a state there are 
no invasions at all. Consequently, the ratio of populations is flxed during the evolution 
meanwhile the system tends towards an uncorrelated spatial distribution for Xg > 0. 
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For all the four models there exist several three-species associations whose members 
invade cyclically each other. For example, in model A one of the the corresponding 

fc) 

states is denoted as $024 where the subscripts refer to the three species being present 
with the same average concentrations, Po = P2 = P4 = 1/3. The corresponding spatial 
distribution is maintained by the cyclic invasions yielding a short average life time for 
the individuals (the survival probability of individuals decreases exponentially with a 
relaxation time r ~ 1.8 MCS) and the correlation length ^ 2.5) is characterizing the 
typical domain size [HI [T31 122] • In general, the stability of the possible three-species 
associations against the external invaders is determined by the food web topology |19j . 
For the model A the food web contains only two equivalent three-species cycles which are 
called defensive alliances because their members protect cyclically each other against the 

(c) 

external invaders. For example, within the state $024 species 2 can only be invaded 
by the "external" species 1. During the cyclic invasion processes both the internal 
species 2 and the substituted invader 1 are killed by their common predator (0), in the 
meantime the third member (4) of the association feeds the species 2 and blocks the 
spreading of invaders. Thus, sooner or later the offspring of invaders will be eliminated 
from the system. One can easily check that other types (3 or 5) of invaders become 
extinct in the same way because of the cyclic symmetry in this food web. For the food 
web B the system has only one cyclic defensive alliance ($135) that will dominate the 
final stationary state although this food web has four additional three-species cycles 
with lower stability. Model C has eight equivalent three-species cycles, and neither of 
them can be considered as a defensive alliance in the sense mentioned above. 

Model D represents another situation where one can find four three-species and 

(c) 

three four-species cycles. In this case the four-species cycle $0234 proved to be the most 
stable spatial association that satisfies the conditions to be a cyclic defensive alliance 

For Xs = 0, the stability of the different states is already determined in a previous 
study when considering the effect of mutations on the population in similar six-species 
predator-prey models fHl- The present analysis will be focused on the effect of site 
exchange between the neutral pairs [X^ > 0) for the absence of mutation. For small 
sizes all the above mentioned states can be observed as a final state. For sufficiently large 
system sizes, however, the system tends towards the real stationary state independent 
of the initial state. 

3. Four-species subsystems 

In this section we study the behavior of several four-species subsystems whose 
understanding will help the interpretation of results obtained for the more complicated 
systems. More precisely, our analysis is focused on three systems (labelled as a, b, and 
c) whose food web are given in Fig. |21 

The food web of subsystem a is equivalent to a four-species cyclic predator-prey 
model where each species has only one predator and one prey. Such a situation can 
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Figure 2. Food webs for the four-species subsystems, 
occur in the models A, B, and D if two species are missing, and its analysis becomes 

(c) 

particularly important for model D where $0234 is the stable state for Xg = 0. In 
addition, similar cyclic dominance has been reported very recently by Traulsen et al. 
|231 121] considering a four-strategy prisoner's dilemma game. 

Some features of this system is already investigated previously by several authors 
for = j^|2niE]- In this case the cyclic invasions maintain a self-organizing pattern 
with Po = pi = P2 = Ps = 1/4 and with a probability P„ = 0.0518(5) to find neutral 
pairs on two neighboring sites. Here P„ characterize the concentration of interfaces 
separating neutral domains (e.g., $9^'' and $2^''). Along these interfaces the evolution is 
blocked until one of the invading predators reaches the interface. For Xs > 0, however, 
along these neutral interfaces the site exchange increases the value of P„ and supports 
the formation of the well-mixed phase of neutral species. 

The MC simulations indicate clearly how the value of Pn increases monotonously 
with Xs until a threshold value [X^ < X^^""^ = 0.02662(2)]. At the same time the 
probability Ppp of predator-prey pairs is decreasing as shown in Fig. IHl 
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Figure 3. Probability of neutral (pluses) and predator-prey (open diamonds) pairs as 
a function of Xg for the four-species cyclic food web (see a in Fig. [5J ■ 

Above this threshold value {Xg > Xj:^"'^) the system segregates into two types ($02^ 
and $S3^) of growing domains which are mixtures of neutral pairs. During the domain 
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growing processes these phases are in contact throughout a wide boundary layer (of type 
$0123) which serves as a species reservoir and sets the compositions to be symmetric |21j . 
Finally, one of the two equivalent mixed states will prevail the whole finite system and 
then Ppp = and P„ = 0.5. 

Subsystem b represents the situation when the three-species cyclic defensive alliance 
($223) is attacked by an external invader (of type 0). The MC simulations indicate 
clearly that after some relaxation process the species dies out and the system remains 
in the state $^23 if < X^^''^ = 0.0527(1). Conversely, the system develops into the 
symmetric mixed phase $02'* (po = P2 = 0.5) for Xs > Xj^^''\ In the vicinity of the 
transition point the visualization of the species distribution illustrates the formation of 
domains of $02^ whose area increases (decreases) slowly above (below) the threshold 
value. The average velocity of the interfaces (separating the domains of $02^ and $123) 
seems to be proportional to Xg — Xj^^''^ as it is observed for another model |2I]. This 
feature can explain why the ordering process becomes so slow in the close vicinity of 
the transition point. 

The four-species subsystem c possesses two equivalent, three-species cycles that 
have two common species (1 and 3). In this case the spatial distribution of species 
segregates into two types ($123 and $103) of growing domains and finally the whole 
system is dominated by one of them with equal probabilities. Notice that the bulk of 
these phases are not influenced by the site exchange mechanism whose effect is limited 
to the boundaries separating the phases ^123 and $103- Furthermore the species within 
the mixed two-species phase $02^ are not able to protect each other because species 1 can 
invade their territories without restraint. This is the reason why the domain growing 
process is not prevented for Xg > and the site exchange causes only a slight variation 
in the velocity of growing process. 

4. Model A 

As mentioned above the six-species model A has two equivalent three-species cycles 
consisting of disjoint sets of species, namely $024 and $135. According to the simulations 
one of these phases (with equal probability) will overwhelm the whole finite system after 
a domain growing process if the site exchange probability does not exceeds a threshold 
value, i.e., Xg < Xj^^"^^ = 0.05592(1). On a sufficiently large scale this domain growing 
process is similar to those that one can observe for the kinetic Ising model below the 
critical temperature if the system is started from a random initial state (for a survey of 
domain growth see ^51). 

Above the threshold value (Xg > X^!^^^) the system develops into one of the 
three (equivalent) two-species states composed from neutral species. The topological 
features of the food web A implies that the phases $03^ ^^'^ ^25 can also be 

considered as defensive alliances because their members protect each other against the 
external predators in the well-mixed state. During the transient time one can observe a 
domain growth with these three types of domains. Apparently this growth is similar to 
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those observed in the three-state Potts model below the critical temperatures l?7j. 
However, along the interfaces separating two domains (say $03^ and ^14) there occur 

(c) 

a four-species cycle $0134 whose transversal extension is limited by invasions from both 
sides. Notice that the corresponding four species subsystem (see the food web a in Fig,|2I) 
exhibit a transition from the four-species cycle ($0134 mentioned example) to one 

of the mentioned neutral pair states at X, = X^J^l As X^^A) > ^^4a) therefore the 
coexistence of the four species is limited to the narrow strips separating the well-mixed 
phases of the neutral pairs. 

The three-color maps have special points (called three-edge vertices) where three 

(c) 

domains meet. Instead of it, here there are patches with phases of either $024 
$^35. This inevitable coexistence of many different phases (during the decomposition 

processes) may be the reason why X(,6^) > X^J^\ 



5. Model B 



In this model a very interesting behavior is found by the MC simulations. One 
can observe four different phases and three subsequent phase transitions when Xg 
is increased. It is more surprising that the transition points are very close to each 
other. To demonstrate the main features Fig. HI illustrates the average values of species 
concentrations in the stationary states reached by the system starting from a random 
initial state for different values of Xg. Here the open (closed) symbols are used for the 
species with odd (even) labels to reduce the confusion coming from the coincidence of 
several data. Furthermore, the neutral pairs are denoted by the same type of symbols. 
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Figure 4. Concentration of species ws. Xg. The MC data are denoted by closed circles, 
open diamonds, closed squares, open circles, closed diamonds, and open squares for 
the species labelled from to 5. The arrows indicate the three transition points. 

The cyclic defensive alliance $^35 is proved to be the most stable solution if 
Xs < X^^^^ = 0.06155(5). In the next region [X^^^^ < X, < X^^^ = 0.06545(5)] 
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the final stationary state will be equivalent to one of the three well-mixed states of 
neutral pairs ($03'' 5 ^ia\ ^25)- the six species survive in the third region of 
Xg in such a way that po = P2 = Pa and pi = ps = ps- If Xg > xj.^^^ = 0.07372(2) 
then only the species with even labels are present and their coexistence with the same 
concentration is maintained by the cyclic invasions (state $024)- 

The transitions at X^i and Xc2 can be considered as a first order transition 
meanwhile the third one exhibits the features of the continuous (critical) transitions. 
The concentration of the species 1, 3, and 5 vanishes continuously when approaching 
Xc3 from below. More precisely, our data are consistent with a power law decrease in 
the close vicinity of the transition point, that is pi = p3 = ps ~ {Xc3 — XsY with 
{3 = 0.37(4). Similar power law behavior (with the same exponent) is found for P„ as 
shown in Fig. El At the same time the MC data show diverging fluctuations for these 
quantities that is another relevant feature of the critical transitions. Unfortunately, we 
could not deduce an adequate (7) exponent to describe this divergency because of the 
large statistical uncertainties in these quantities. 

10-' r 




Figure 5. Log- log plot of the MC data (open squares) for the probability of neutral 
pairs vs. -^^^3^'' ~ -^s- The fitted power law is denoted by the straight line with a slope 
of 0.373. 

In the presence of absorbing state (s) the extinction of a single species in the spatial 
models has very robust features and the corresponding transition belongs to the directed 
percolation (DP) universality class characterized by a higher exponent {Pdp = 0.58(2) 
for the two-dimensional systems) |28|, 123 IHOl 1^ ■ In the present model, however, three 
species die out simultaneously and the background can mediate some interactions among 
the vanishing species. Evidently, further analysis is necessary to quantify the features 
of this transition as well as to clarify the reasons causing the distinct behavior. 

In contrary to the model A here the formation of one the three equivalent well- 
mixed phases of neutral species is limited to a narrow region [X^l^^ < Xs < X^2^^ = 
0.06545(5)]. The formation of these structures differs significantly from those described 
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for model A because the food web B does not provide mutual protection for a neutral pair 
against another neutral pair. Instead, the well-mixed phases of the neutral pairs destroy 
cyclically each other ($03'' destroys $^4^ destroys $25^ destroys $03'*) ^ complicated 
process. For example, the offspring of a single species can occupy the whole territory 
of the phase $^4 . As a result, if the growing domains of the phases $03^ and $^4 meet 
along a common boundary then the domain $^4 is transformed (via the invasion) into 
a domain of $0'^'' that is also unstable against the invasions of species 2 and 5. 

Within the second region of Xg the six species can coexist for a sufficiently long time 
in such a structure where the species with even labels are in minority. The formation 
of the mixed phases of the neutral species begins in those patches where two species 
of minorities are missing due to the fluctuations. According to the previous analysis of 
the four-species subsystem b the spontaneous formation of the well-mixed phase of the 
corresponding neutral species is permitted because the condition Xg > Xj^^^'' is fulfilled 
within the given region of Xg. By this means sufficiently large domains of the states 
$03^ $14^ and $25'' happen to form by chance. During their expansion these domains 
meet one of their rivals and then one of them is destroyed as described above. 




Figure 6. Time dependence of tlie probability of neutral (dashed line) and predator- 
prey (solid line) pairs for Xg = 0.0654 and L — 2000. The MC data are smoothed out 
by averaging over 2000 MCS. 

Figure ini shows the variation of P„(t) and Ppp{t) for a given Xg close to the second 
transition point. At the beginning these curves are smooth. Later, however, one can 
observe larger and larger variations reflecting the destruction of domains with types 
$03^ $14^ and $25^ whose average size increases monotonously. The series of these 
larger and larger invasion processes is ended when only one of these phases prevail the 
whole system. 

It is conjectured that this (rock-scissors-paper like) cyclic dominance between the 
two-species associations plays a crucial role in the maintenance of the six-species phase. 
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6. Model C 

This model has eight equivalent three-species cyclic associations as mentioned above. 
The simulations for a system size of L = 400 have justified that any of them can be the 
final stationary state independently of the value of X^. 

It is expected that this "mono-domain" state is developed via a domain growing 
process. This phenomenon can be described by the variation of the probability of neutral 
pairs vanishing in the bulk of "mono-domain" states. The results of MC simulations 
(for L = 2500) are plotted in Fig. [7| where the data are averaged over a sampling time 
increased from 2 to 2000 MCS. 
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Figure 7. Log-log plot of the probability of neutral pairs vs. time for Xs = 0.05 
(upper solid line) and Xg — (central solid line). The lower solid line shows the same 
quantity in the model A for Xg — 0- To compare the MC results with the theoretical 
expectations, the upper dashed line represents the function P„(t) = 0.12/ ln(t) while 
the lower one shows an algebraic decay, P„(t) = 0.25/v^. 

Figure [7| shows that the very slow domain growing process for Xg = becomes even 
slower for Xg = 0.05. Indeed we cannot exclude the blocking of coarsening process as it 
happens in the voter model for higher dimensions {d > 2) Keeping in mind the large 
statistical uncertainties (for large times) these data remind us to the growing processes 
observed for the two-dimensional voter model where the concentration of domain walls 
vanishes as l/ln(t) plotted also in Fig. [3 |38| l82] . The voter model exemplifies the 
coarsening process in a broad class of models undergoing a phase ordering without 
surface tension |23j . In the present model the definition of the boundary separating two 
cyclic three-species phases is confusing because these phases may involve one or two 
common species (compare ^qi^ with ^qi2 or $123)- 

Strikingly different growth process is found for model A (discussed previously) 
where the system has only two cyclic three-species phases composed from disjoint sets 



Competing associations in six-species predator-prey models 12 

of species ($024 ^135)- order to visualize the relevant differences Fig. shows also 
the MC results obtained in the model A for Xg = and L = 4000. These numerical 
results indicate that here the asymptotic behavior of Pnit) is close to the theoretical 
prediction {Pn{t) ~ l/v^) characterizing the growth controlled by the reduction of an 
interfacial energy 12Z1 • 

Evidently, further research is required to clarify those ingredients of these models 
that are responsible for the different growth processes. Furthermore, here it is worth 
mentioning that the logarithmic growth supports the maintenance of biodiversity for a 
very long time if the system is sufficiently large. 



7. Model D 

For this food web the system has a four-species cychc defensive alliance ($0234 &s shown 
in Fig. ^ that is stable if = ^Hj- According to the simulations for weak rate of site 
exchange the species 1 and 5 die out within a short transient time. In the corresponding 
four-species subsystem (see Sect. Ej) the cyclic invasions can sustain their coexistence if 
Xs < Xj^^"-\ For higher Xg the spatial distribution of these four species decomposes into 
two types of two-species domains ($03'' '^24'') &s described above. If Xg is increased 
then this process becomes faster, and above a threshold value the domains of phases 
$03'' and $24'' built up by these processes before their mortal enemy (species 1 and 
5) would die out. One can easily check that the species 1 (5) can occupy the whole 
territory of the phases $03'' {^24)- Thus, for a sufficiently intensive mixing the species 
1 and 5 will have an enhanced chance to survive because they are fed by the emerging 
domains of $03'' and $24'*- large Xg this mechanism can be so effective that finally 
the phase $^5'' will dominate the whole system as shown in Fig. |H| Notice that the well 
mixed phase $15^ is also a defensive alliance because species 1 and 5 protect each other 
against the different invaders. 

Figure |S1 shows that four types of final stationary states can be distinguished when 
Xg is increased. In the first region [Xg < X^.^^^ = Xjij:"-^) the species of the phase $0234 
take place with the same concentration meanwhile P„ and Ppp vary with Xg as plotted in 
Fig. El Starting from a random initial state this system evolves into one of the neutral- 
pair phases $[,3^ and $24'* if ^d^^ < Xg < X^^^^^ where the value of X^^^^^ depends on 
the size L as it will be discussed below. Within the third region {X^^^^^ < Xg < X^^^^) 
the six species coexist in such a way that po = P2 = Ps = P4 > Pi = P5 > 0. If 
Xg > X^^^^ = 0.0581(1) then only the species 1 and 5 will survive with forming a phase 

^15 • 

In this model the first and the third transitions exhibit a sudden change both in 
concentrations and in pair probabilities (P„ and Ppp)- However, the classification of the 
second transition (as well as the numerical value of X^2^^) is difficult because of its 
unusual features. Namely, pi and ps decrease monotonously when Xg is decreased 
within the third region. At the same time we have observed a dramatic increase 
in the fiuctuation of concentration for the majority species (see Fig. El meanwhile 
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Figure 8. Concentration of species vs. Xs in the model D. The MC data are denoted 
by symbols with fourfold symmetry (pluses, Xs, open squares, and diamonds) for those 
species belonging to four-species cyclic defensive alliance, and by open circles for the 
species 1 and 5. The arrows indicates the critical values oi Xs where transitions appear. 



the corresponding average concentrations are close to 1/4 (see Fig. IH}. To avoid the 
undesired effects of these fluctuations the simulations were performed for such a large 
system sizes (the linear size is increased up to L = 2600) where the current value of 
concentrations were significantly larger than the amplitude of fluctuation characterized 
by their standard deviation dependent on size [^J~xJn). For smaller sizes the system 

evolves into the phase $03'' or $24^- 
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Figure 9. Fluctuation of concentration as a function of Xg for the species as denoted 
by the figures. 



Figure ini demonstrates clearly the absence of fluctuations for the missing species as 
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well as for the phases composed of neutral species. The MC data refer to an exponential 
increase in xo: X2, Xs? X4 if approaches Xc2 from above. 

Within the region of large fluctuations the time dependence of concentrations shows 
the occurrence and destruction of domains with the structure $03^ and $24^- These 
domains are invaded by the sparsely dispersed species 1 and 5. The resultant phase, 
however, is not stable at the given Xg and sooner or later the phases $03'* or $24'* 
formed again. In some sense the evolution of the spatio-temporal pattern is similar to 
those described by the mentioned spatial rock-scissors-paper games [HI [ini EI as well 
as in the forest-fire models ISHl EZI- In the present case two cycles interfere in a 
complex way because cyclic dominance emerges between the different set of associations 
and this mechanism supports the maintenance of biodiversity with more associations 
(and species). 

Due to the mentioned difficulties we couldn't determine accurately how 
concentrations pi and ps vanish when approaching X^2^'* from above. In the light 
of the MC results one can conclude a weakly first order phase transition (sudden jump 
in pi = Ps) accompanied with an increase in xi = Xb illustrated in Fig. El 

8. Conclusions 

The numerical analysis of the four different six-species, spatial predator-prey models 
has highlighted the large variety of states and phenomena (transitions) that can occur 
in self-sustaining, multi-species ecological systems. Due to their simplicity the present 
models have allowed us to consider quantitatively some interesting phenomena. Instead 
of recalling the curious features found for these oversimplified models henceforth I wish 
to concentrate on those general messages concluded from these examples. 

One of the relevant messages is related to the increase in the degree of freedom. 
Many recent works are focused on spatial systems with one, two, or three species. For 
larger number of species the present models exemplify how the simple dynamical rules 
build up more complex objects (alliances) that can be considered as additional species 
participating in the evolutionary competition. Shortly, the degree of freedom (here 
the number of species) is increased by the spontaneously occurring alliances having 
particular composition and spatial structure. One can believe that the food web of 
these systems should be extended by these complex objects (with a suitable invasion 
rates) to have a better picture about the whole spatial system. 

In the context of game theory (for details see the refs. [H El El IH IE] ) , the interaction 
of species is described by the payoff matrix (the adjacency matrix of the food web in 
the present models) that determines the Nash equilibrium from which the unilateral 
deviation yields loss. For many evolutionary (dynamical) rules the Nash equilibrium 
represents an evolutionary stable state under the mean-field conditions. The Nash 
equilibrium is not necessary unique and it can be a mixed state (strategy). In the 
spatial models with short range interactions several species can be absent from a given 
territory and the Nash equilibrium for the corresponding subsystem can be distinct from 
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those characterizing the complete system. The noises affect the stabihty of these spatial 
states |HH1 and modify the average velocity of invasions along the interfaces separating 
the competing phases [21] • The above results justify that the invasion processes between 
the "stable" states of subsystems can play a crucial role in the formation of the final 
stationary states. 

For many cases the final stationary state is reached through a coarsening 
process where the alliances (degenerated Nash equilibria) form growing domains. 
The corresponding domain growth process is affected by the additional phases (and 
mechanisms) along interfaces. It is not yet known how these effects modify the growth 
process yielding significantly slower (logarithmical) domain growth. 

Most of the (first order) transitions occur at such a value of Xg where the invasion 
rate between two dominant alliances vanishes as detailed previously [21]. Remarkable 
exceptions are the fourth transition in model B and the second one in model D. In both 
cases more than one species die out simultaneously and their extinction is accompanied 
by an enhanced concentration fluctuation for the surviving species and the resultant 
spatio-temporal structures can mediate unusual interactions between the vanishing 
species. These transitions are distinct inherently and their classification requires further 
time-consuming numerical investigations. 

The spontaneous formation of complex spatio-temporal structures is one of the 
most interesting predictions of these simple models because it sustains an abundance of 
species. Furthermore, this feature mixes the concept of the part and whole. At the same 
time the corresponding system can be interpreted as a set of objects that create objects 
in a way characterizing the living systems [SH] • In the light of the present results one can 
easily imagine that the choice of more complex food webs can result in the emergence 
of more and more complex objects (e.g., alliances of alliances with particular spatio- 
temporal structures). Besides it these models imply the possibility that the evolution 
of food web by mutations [101 HH 112] and the evolution of the spatial ecological system 
are strongly related to each other while the spatial inhomogeneity can be enhanced in 
the system. 
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